% code used to create Fig. 7 in "An alternative pathway for delivery of 
% power by outer hair cells to cochlear traveling waves", Samaras and Meaud
% evaluate power transfer to the BM
% last modified 10/30/2025
% Julien Meaud
% julien.meaud@me.gatech.edu
% ------------------------------------------------------------------------

close all 
clear all

computeFluidForceFromGlobalMatrix=1;
calculateUsingGlobalMatrices=0; % set to 1 if you need to calculate power and force (when =0, these are already in datafile)

plotIndB=0;

CF=20000; % CF of the location used in the plots

setPlottingOptions

modelV=[1:4];

iModelPlot=0;
for iModel=modelV

    iModel    
    iModelPlot=iModelPlot+1;
       
    ls='-';
    lw=2;

    % load the data ----------------------------------------------
    if iModel==1
        baseFileName='PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25'; 
    elseif iModel==2
        baseFileName='PTLFD_G35_Act3.15B1_15A_rEps3_1_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25'; 
     elseif  iModel==3
        baseFileName='PTLFD_G35_Act1.65B1_15A_rEps3_1_rKdcKohc_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25';   
    elseif iModel==4
        baseFileName='PTLFD_G35_Act3.0B1_15A_rEps3_1_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25';            
    end
        
    if calculateUsingGlobalMatrices==1            
        global L NmodesBM NmodesTM
        global nMechTotal nElecTotal
        global Nme CHANGE_ME_1DOF xMesh  randomNumber_TMmass randomNumber_TMarea % to generate global Ks matrix
        global Nrw Nd 
        global elecLongCoupling  INCLUDE_ST_CABLES INCLUDE_SM_CABLES INCLUDE_SV_CABLES nHelicTotal
        
        computeFluidForceFromGlobalMatrix=1;     
        plotIndB=0;        
        
        NmodesBM=1;
        NmodesTM=2;
        
        elecLongCoupling =1;
        INCLUDE_ST_CABLES=1;
        INCLUDE_SM_CABLES=1;
        INCLUDE_SV_CABLES=1;
        fileName=['DataWithFluidPressure/',baseFileName,'_withPressure.mat'];

    else
        fileName=['PowerAnalysis/',baseFileName,'_withForcePower.mat'];
    end

    data=load(fileName);
    % ---------------------------------------------------------------------

    titleStr=['M',int2str(iModel)];
    iPlot=iModel;



    [~,indexFreq] =min(abs(CF-data.StimulusOptionsStruct.freq));        
    
    ActLevel=data.ModelOptionsStruct.electrical.actLevel;

    xMesh=data.ModelGeneratedVarsStruct.xMesh;
    freq=data.StimulusOptionsStruct.freq;

    [~,indexPeak]=max(abs(data.ModelOutputsStruct.structural.uBM(:,indexFreq)));
    
    uBM=data.ModelOutputsStruct.structural.uBM(indexPeak,:);
    uTMR=data.ModelOutputsStruct.structural.uTMR(indexPeak,:);
    uTMB=data.ModelOutputsStruct.structural.uTMB(indexPeak,:);
    uDC=data.ModelOutputsStruct.structural.uDC(indexPeak,:);

    phiSV=data.ModelOutputsStruct.electrical.phiSV(indexPeak,:);
    phiSM=data.ModelOutputsStruct.electrical.phiSM(indexPeak,:);
    phiOHC=data.ModelOutputsStruct.electrical.phiOHC(indexPeak,:);
    phiST=data.ModelOutputsStruct.electrical.phiST(indexPeak,:);

    vBM=uBM.*1i.*freq*2*pi;
    vTMR=uTMR.*1i.*freq*2*pi;
    vTMB=uTMB.*1i.*freq*2*pi;
    if data.ModelOptionsStruct.structural.BMlongCoupling==1
       uBMrot=data.ModelOutputsStruct.structural.uBMrot(indexPeak,:);
        vBMrot=uBMrot.*1i.*freq*2*pi;
    end

    x=xMesh(indexPeak);
    uLoc=[uBM;uTMR;uTMB;uDC];
    phiLoc=[phiSV;phiSM;phiOHC;phiST];
    
      
       
   if calculateUsingGlobalMatrices==1
       fileNameOutput=[];
        [structForces,structPower]=forceAnalysis(data,fileNameOutput);

       
        ModelGeneratedVarsStruct=data.ModelGeneratedVarsStruct;
        ModelOptionsStruct=data.ModelOptionsStruct;
        ModelOutputsStruct=data.ModelOutputsStruct;
        ModelOutputsStruct.pressure=rmfield(ModelOutputsStruct.pressure,'pFull');
        ModelOutputsStruct.pressure=rmfield(ModelOutputsStruct.pressure,'P_box_SpecifiedFreqs');
        OptionalStruct=data.OptionalStruct;
        SaveOptionsStruct=data.SaveOptionsStruct;
        SimulationMetricsStruct=data.SimulationMetricsStruct;
        StabilityAnalysisStruct=data.StabilityAnalysisStruct;
        StimulusOptionsStruct=data.StimulusOptionsStruct;

         fileName=['PowerAnalysis/',baseFileName,'_withForcePower.mat'];
         save(fileName,'ModelGeneratedVarsStruct','ModelOptionsStruct','ModelOutputsStruct','OptionalStruct','SaveOptionsStruct','SimulationMetricsStruct',...
             'StabilityAnalysisStruct','StimulusOptionsStruct','structForces','structPower') 
   else
       structForces=data.structForces;
       structPower=data.structPower;
   end

    FbmPeak=structForces.FbmMat(indexPeak,:);
    FbmViscousPeak=structForces.FbmViscousMat(indexPeak,:);
    FfluidPeak=structForces.FfluidMat(indexPeak,:);
    FdcPeak=structForces.FdcMat(indexPeak,:);
    FpcPeak=structForces.FpcMat(indexPeak,:);

    Power_bm_on_BM_global=structPower.Power_bm_on_BM_global(indexPeak,:);
    Power_bmViscous_on_BM_global=structPower.Power_bmViscous_on_BM_global(indexPeak,:);
    Power_dc_on_BM_global=structPower.Power_dc_on_BM_global(indexPeak,:);
    Power_pc_on_BM_global=structPower.Power_pc_on_BM_global(indexPeak,:);
    Power_fluid_on_BM_global=structPower.Power_fluid_on_BM_global(indexPeak,:);

    Lp=20;
    pref=20e-6*10^(Lp/20);
    
    figure(1)
    subplot(2,length(modelV),iModelPlot)
    plot(freq/1000,20*log10(abs(uBM*1e7/pref)),'Color','b','LineStyle',ls,'LineWidth',lw)
    hold on
    plot(freq/1000,20*log10(abs(uDC*1e7/pref)),'Color','r','LineStyle','--','LineWidth',lw)
    ylabel('Disp. re P_{ec} (dB 1 nm/Pa)')
    if iPlot==1
        legend('BM','DC')
    end
     xlim([3 30])
    ylim([0 100])
    title(titleStr)

    subplot(2,length(modelV),length(modelV)+iModelPlot)
    plot(freq/1000,(angle(uDC./uBM))/(2*pi),'Color','r','LineStyle',ls,'LineWidth',lw)
    hold on
    ylabel('Phase u_{DC} re u_{BM} (cycles)')
    legend('Original','\alpha_{C_{dc}}=0','\alpha_{C_{dc}}=0.5')

    ylim([-0.5 0.5])
    title(titleStr)
    xlim([3 30])
     
       
        
          
       fig56=figure(56);
        set(fig56,'units','inches','Position',[2 2 7 6])  
        subplot(4,length(modelV),iModelPlot)
        FrefGlobal=max(abs(FbmViscousPeak(50:end)));            
        if plotIndB==1
            plot(freq/CF,20*log10(abs(FfluidPeak/FrefGlobal)),'Color','k','LineWidth',2,'DisplayName','F_{fluid/bm}')
            hold on
            plot(freq/CF,20*log10(abs(FdcPeak/FrefGlobal)),'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','F_{dc/bm}')
            plot(freq/CF,20*log10(abs(FpcPeak/FrefGlobal)),'Color','r','LineStyle','--','LineWidth',lw,'DisplayName','F_{pc/bm}')
            plot(freq/CF,20*log10(abs(FbmViscousPeak/FrefGlobal)),'Color','b','LineStyle',':','LineWidth',2,'DisplayName','F_{bm}^v')
            plot(freq/CF,20*log10(abs(FbmPeak/FrefGlobal)),'Color','b','LineStyle','--','LineWidth',2,'DisplayName','F_{bm}^{lc}')
        
            if iModel==1
                ylabel(['Norm. force',char(10),'(dB)'])      
                legend('orientation','horizontal','AutoUpdate','off')
            end
            line([1 1],[-70 10],'Color','k','LineStyle',':')
            ylim([-70 10])
        else
            semilogy(freq/CF,(abs(FfluidPeak/FrefGlobal)),'Color','k','LineWidth',2,'DisplayName','F_{fluid/bm}')
            hold on
            semilogy(freq/CF,(abs(FdcPeak/FrefGlobal)),'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','F_{dc/bm}')
            semilogy(freq/CF,(abs(FpcPeak/FrefGlobal)),'Color','r','LineStyle','--','LineWidth',lw,'DisplayName','F_{pc/bm}')

            semilogy(freq/CF,(abs(FbmViscousPeak/FrefGlobal)),'Color','b','LineStyle',':','LineWidth',2,'DisplayName','F_{bm}^v')
            semilogy(freq/CF,(abs(FbmPeak/FrefGlobal)),'Color','b','LineStyle','--','LineWidth',2,'DisplayName','F_{bm}^{lc}')
        
            if iModel==1
                ylabel(['Norm. force'])      
                legend('orientation','horizontal','AutoUpdate','off')
            end
            line([1 1],[0.0003 3],'Color','k','LineStyle',':')
            ylim([0.003 30])
        end
            
            xlim([0.3 1.5])            
            title(titleStr)
            set(gca,'XTick',[0.5 1 1.5],'XTickLabel',[])

            subplot(4,length(modelV),length(modelV)+iModelPlot)
            patch(xPatch,yPatch1,orange,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])
            hold on
            patch(xPatch,yPatch2,green,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])
            patch(xPatch,yPatch1bis,orange,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])            
            plot(freq/CF,unwrap(angle(FfluidPeak./vBM))/(2*pi),'Color','k','LineWidth',2,'DisplayName','fluid')
            hold on
            plot(freq/CF,unwrap(angle(FdcPeak./vBM))/(2*pi),'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','DC pathway')
            plot(freq/CF,unwrap(angle(FpcPeak./vBM))/(2*pi),'Color','r','LineStyle','--','LineWidth',lw,'DisplayName','PC pathway')
            plot(freq/CF,unwrap(angle(FbmViscousPeak./vBM))/(2*pi),'Color','b','LineStyle',':','LineWidth',2,'DisplayName','F_{bm}^v')
            plot(freq/CF,unwrap(angle(FbmPeak/FrefGlobal./vBM))/(2*pi),'Color','b','LineStyle','--','LineWidth',2,'DisplayName','F_{bm}^{lc}')
            
            ylim([-0.51 0.51])
            if iModel==1
                ylabel(['Phase F re v_{BM}',char(10),'(cycles)'])
            end  
            xlim([0.3 1.5])
            %xlabel('Norm. frequency')
            set(gca,'XTick',[0.5 1 1.5],'XTickLabel',[])
            line([1 1],[-0.5 0.5],'Color','k','LineStyle',':')

            powerRef=max(abs(Power_bmViscous_on_BM_global));%+Power_bm_on_BM_global));
            subplot(4,length(modelV),2*length(modelV)+iModelPlot)
            plot(freq/CF,Power_fluid_on_BM_global./powerRef,'Color','k','LineWidth',2,'DisplayName','\Pi_{fluid/bm}','LineWidth',2)
            hold on
            plot(freq/CF, Power_dc_on_BM_global./powerRef,'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','\Pi_{dc/bm}')
            plot(freq/CF, Power_pc_on_BM_global./powerRef,'Color','r','LineStyle','--','LineWidth',lw,'DisplayName','\Pi_{pc/bm}')
            plot(freq/CF, (Power_bmViscous_on_BM_global)./powerRef,'Color','b','LineStyle',':','LineWidth',2,'DisplayName','\Pi_{bm}^v')      
            hold on
            plot(freq/CF, (Power_bm_on_BM_global)./powerRef,'Color','b','LineStyle','--','LineWidth',2,'DisplayName','\Pi_{bm}^{lc}')     
                       
            Power_tot=Power_bm_on_BM_global+Power_bmViscous_on_BM_global+Power_fluid_on_BM_global+ Power_dc_on_BM_global+ Power_pc_on_BM_global;
           % plot(freq/CF, (Power_tot)./powerRef,'Color','k','LineStyle','--','LineWidth',2,'DisplayName','\Pi_{bm}^{tot}')     
            xlim([0.3 1.5])
            
             if iModel==1
             
                ylabel(['Norm. power',char(10),'delivered to BM'])
                legend('orientation','horizontal','AutoUpdate','off')
             end
         
           set(gca,'XTick',[0.5 1 1.5],'XTickLabel',[])
           line([1 1],[-15 20],'Color','k','LineStyle',':')         
            
            subplot(4,length(modelV),3*length(modelV)+iModelPlot)
            Zbm_global=FfluidPeak./vBM;

            ZbmLC_global=FbmPeak./vBM;
            plot(freq/CF,real(Zbm_global)*1e-1*1e3,'Color','k','LineWidth',2,'DisplayName','Fluid')
            hold on
            xlim([0.3 1.5])
            line([freq(1),freq(end)]/CF,[0 0],'Color','k','LineStyle',':')
            xlabel('Freq. re BF')
            if iModel==1
                ylabel(['BM resistance',char(10),'mN/s/m'])
            
           end
           set(gca,'XTick',[0.5 1 1.5],'XTickLabel',[0.5 1 1.5])
           line([1 1],[-1 1]*10,'Color','k','LineStyle',':')

           fig57=figure(57)
           subplot(2,4,iModel)
           
            plot(freq/CF, Power_dc_on_BM_global./powerRef,'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','\Pi_{dc/bm}')
            hold on
            plot(freq/CF, -Power_pc_on_BM_global./powerRef,'Color','r','LineStyle','--','LineWidth',lw,'DisplayName','-\Pi_{pc/bm}')
            legend
            xlim([0.3 1.5])  
            line([1 1],[-15 20],'Color','k','LineStyle',':')
            
            subplot(2,4,4+iModel)       
            plot(freq/CF, Power_pc_on_BM_global./Power_dc_on_BM_global,'Color','m','LineStyle',ls,'LineWidth',lw,'DisplayName','\Pi_{dc/bm}')
            hold on  
            legend
            xlim([0.3 1.5])
            line([1 1],[-15 20],'Color','k','LineStyle',':')
end

width=0.19;
height=0.19;
leftX=0.09;
gapX=0.03;
nRow=4;
figure(56)
for iModel=modelV
    if length(modelV)==4
        set(subplot(nRow,length(modelV),iModel),'Position',[leftX+(width+gapX)*(iModel-1) 0.77 width height]);
    else
        subplot(nRow,length(modelV),iModel)
    end
    if iModel>1
        if plotIndB==1
            set(gca,'YTick',[-80:20:0],'YTickLabel',[])
        else
            set(gca,'YTick',[0.001 0.01 0.1 1 10],'YTickLabel',[])
        end
    else
        if plotIndB==0
        set(gca,'YTick',[0.001 0.01 0.1 1 10])
        end
    end
     if plotIndB==1
        text(0.35,-0,char('A'+iModel-1),'FontWeight','bold','FontSize',16)
     else
         text(0.35,10,char('A'+iModel-1),'FontWeight','bold','FontSize',16)
     end
end


for iModel=modelV
    if length(modelV)==4
        set(subplot(nRow,length(modelV),length(modelV)+iModel),'Position',[leftX+(width+gapX)*(iModel-1) 0.53 width height]);
    else
        subplot(nRow,length(modelV),length(modelV)+iModel)
    end
    if iModel>1
        set(gca,'YTick',[-0.5 -0.25 0 0.25 0.5],'YTickLabel',[])
    else
        set(gca,'YTick',[-0.5 -0.25 0 0.25 0.5])
    end
    text(0.35,0.4,char('E'+iModel-1),'FontWeight','bold','FontSize',16)

    if iModel==1
        t1 = text(0.32,0.05,'Power added','Color','g','FontSize',8);
        t2 = text(1,0.4,['Power',char(10),'removed'],'Color','r','FontSize',8);
        t3 = text(0.4,-0.4,'Power removed','Color','r','FontSize',8);
    end
end

for iModel=modelV
    if length(modelV)==4
        set(subplot(nRow,length(modelV),2*length(modelV)+iModel),'Position',[leftX+(width+gapX)*(iModel-1) 0.32 width height]);
    else
        subplot(nRow,length(modelV),2*length(modelV)+iModel)
    end    
 
    patch(xPatch,[-8 -8 0 0],orange,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])
    hold on
    patch(xPatch,[0 0 8 8],green,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])    

    if iModel<=3
        set(gca,'YTick',[-2 -1 0 1 2])
    else
        set(gca,'YTick',[-4 -2 0 2 4])
    end

    if or(iModel==1,iModel==2)
        ylim(([-1.3 2.8]))
        text(0.35,2.4,char('I'+iModel-1),'FontWeight','bold','FontSize',16)
    
    elseif iModel==3
        ylim([-2.8 3.6])
        text(0.35,3,char('I'+iModel-1),'FontWeight','bold','FontSize',16)
    elseif iModel==4
        ylim([-5.2  5.4])
         text(0.35,4,char('I'+iModel-1),'FontWeight','bold','FontSize',16)
    end
 
    if iModel==1
        t4 = text(0.35,0.9,['Power',char(10),'added'],'Color','g','FontSize',8);
        t5 = text(0.35,-0.5,['Power',char(10),'removed'],'Color','r','FontSize',8);
    end

end

for iModel=modelV
    if length(modelV)==4
        set(subplot(nRow,length(modelV),3*length(modelV)+iModel),'Position',[leftX+(width+gapX)*(iModel-1) 0.07 width height]);
    else
        subplot(nRow,length(modelV),3*length(modelV)+iModel);
    end
   
    if iModel==2
        set(gca,'YTick',[-7.5 -5 -2.5 0 2.5])
        ylim([-8 2.5])
        text(0.35,-6.2,char('M'+iModel-1),'FontWeight','bold','FontSize',16)        
        hold on
        
    elseif iModel==3
        set(gca,'YTick',[-2.0 -1 -0.5 0 0.5 1 2])
        ylim([-1 0.7])
        text(0.35,-0.7,char('M'+iModel-1),'FontWeight','bold','FontSize',16)
    else
        set(gca,'YTick',[-2.0 -1 0 1 2])
         ylim([-1.5 1])
          text(0.35,-1.1,char('M'+iModel-1),'FontWeight','bold','FontSize',16)
    end
    xlim([0.3 1.5])

    patch(xPatch,yPatch1_Zbm,orange,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])
    patch(xPatch,yPatch2_Zbm,green,'FaceAlpha',transparency_Patch,'HandleVisibility','off','EdgeAlpha',[0.2])

    if iModel==1
        t1_Zbm = text(1.1,-0.5,['Power',char(10),'added'],'Color','g','FontSize',8);
        t2_Zbm = text(0.55,0.5,['Power',char(10),'removed'],'Color','r','FontSize',8);
    end
    
end

